;+
; NAME:
;       PLFIT
;
; PURPOSE:
;       PLFIT fits a power-law distributional model to data.  Source:
;       http://www.santafe.edu/~aaronc/powerlaws/ 
;       
;       This is an IDL implementation of plfit.m found at the site
;       above.
;
; CATEGORY:
;       Data-Fitting
;
; CALLING SEQUENCE:
;       alpha = PLFIT(x, XMIN=xmin, L=L)
;
; INPUTS:
;       X -- The data to be fit.
;
; OPTIONAL INPUTS:
;       RANGE  -- 
;       SAMPLE -- 
;       LIMIT  -- Hard upper limit on XMIN
;       XMIN   -- Fixed value of xmin
;
; KEYWORD PARAMETERS:
;       FINITE  -- 
;       NOWARN  -- Do not print warning messages.  Error messages will
;                  still be displayed. 
;       NOSMALL -- 
;       LOGDEL  -- For large data sets, speed-up in run time by
;                  searching for best XMIN fit by incrementing in
;                  Delta(Log X) space rather than by individual data
;                  points.  This option is automatically set for data
;                  sets w/ N >= 10^4.  Setting LOGDEL=0 will override
;                  this automatic setting.
;
; OUTPUTS:
;       ALPHA -- Maximum-likelihood power law slope for the input
;                data.
;
; OPTIONAL OUTPUTS:
;       XMIN -- Minimum data value for which the power-law slope is
;               valid, as estimated from the data.
;       L    -- Log Likelihood of the fit.
;
; COMMON BLOCKS:
;       NONE
;
; PROCEDURE:
;       PLFIT(x) estimates x_min and alpha according to the
;       goodness-of-fit based method described in Clauset, Shalizi,
;       Newman (2009). x is a  vector of observations of some quantity
;       to which we wish to fit the power-law distribution p(x) ~
;       x^-alpha for x >= xmin.  PLFIT automatically detects whether x
;       is composed of real or integer values, and applies the
;       appropriate method. For discrete data, if min(x) > 1000, PLFIT
;       uses the continuous approximation, which is a reliable in this
;       regime.
;
;       The fitting procedure works as follows:
;       1) For each possible choice of x_min, we estimate alpha via
;          the method of maximum likelihood, and calculate the
;          Kolmogorov-Smirnov goodness-of-fit statistic D.
;       2) We then select as our estimate of x_min, the value that
;          gives the minimum value D over all values of x_min.
;
;       Note that this procedure gives no estimate of the uncertainty
;       of the fitted parameters, nor of the validity of the fit. 
;
;       Example:
;          x = (1-rand(10000,1)).^(-1/(2.5-1));
;          alpha = plfit(x, XMIN=xmin, L=L)
;
;       The output 'alpha' is the maximum likelihood estimate of the
;       scaling exponent, 'xmin' is the estimate of the lower bound of
;       the power-law behavior, and L is the log-likelihood of the
;       data x>=xmin under the fitted power law.
;
; LICENSE INFORMATION:
;       The plfit.m from which this IDL routine was derived is
;       Copyright (C) 2008-2012 Aaron Clauset (Santa Fe Institute).
;       Distributed under GPL 2.0
;       (http://www.gnu.org/copyleft/gpl.html)
;       PLFIT comes with ABSOLUTELY NO WARRANTY
;
;       This IDL version comes with even less warranty.
;
; MODIFICATION HISTORY:
;
;       Created:  03/23/13, TPEB -- Initial IDL version, derived from
;                                   plfit.m V1.0.11 (2012 January)
;                                   found at http://tuvalu.santafe.edu/ 
;                                   ~aaronc/powerlaws/plfit.m
;       Modified: 04/01/13, TPEB -- Changed method of determining
;                                   whether input is INTS or REAL
;                                   because earlier method could chew
;                                   through some serious memory.
;       Modified: 04/02/13, TPEB -- Removed extraneous print
;                                   statements.
;       Modified: 04/04/13, TPEB -- Updated Documentation.
;       Modified: 04/11/13, TPEB -- Added support for specifing a
;                                   maximum value of x, above which
;                                   the power-law description does not
;                                   apply.  A correction term to dL/da
;                                   is added, mandating numerical
;                                   root-finding for \hat{\alpha}.
;       Modified: 04/16/13, TPEB -- Working to correct above.
;       Modified: 05/07/13, TPEB -- Added /NaN keyword to INTERPOL
;                                   commands to conform with IDL 8
;                                   usage.
;       Modified: 09/03/13, TPEB -- Added LOGDEL keyword to speed
;                                   computation for very large data
;                                   sets by incrementing XMIN in set
;                                   Delta(Log X) steps rather than by
;                                   data point.
;
;-

FUNCTION PLFIT, x, XMIN=xmin, ELL=L, RANGE=vec, SAMPLE=sample, LIMIT=limit, $
                FINITE=finite, NOWARN=nowarn, NOSMALL=nosmall, XMAX=xmax, $
                LOGDEL=logdel, VERBOSE=verbose
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  ;; Set limit where LOGDEL kicks in automatically
  ldmax = 1e4
  
  ;; Parse keywords
  finite  = KEYWORD_SET(finite)
  nowarn  = KEYWORD_SET(nowarn)
  nosmall = KEYWORD_SET(nosmall)
  verbose = KEYWORD_SET(verbose)
  
  IF n_elements(xmin) THEN xminx = xmin
  
  ;; Parse LOGDEL - Y for large sets, but settable to 0b
  logdel = (KEYWORD_SET(logdel) || (n_elements(x) GE ldmax)) $
           && ~(n_elements(logdel) && ~KEYWORD_SET(logdel))
  
  ;; Check IDL version for INTERPOL usage
  oldver = fix(!version.release[0]) LE 7
  
  ;; Error check inputs
  nx = n_elements(x)
  IF nx LT 2 THEN BEGIN
     message,'Error: Input X must be present and contain at least two '+$
             'elements, else I cannot compute a power-law, now can I?  '+$
             'Returning.',/cont
     RETURN,0.
  ENDIF
  
  IF n_elements(xmax) THEN BEGIN
     xlim = 1b
     xind = where(x LE xmax, nlim)
     IF ~nlim THEN BEGIN
        message,'Error: Input XMAX must be larger than the smallest value '+$
                'of X.  Returning.',/cont
        RETURN,0.
     ENDIF
     retx = x
     x = x[xind]
     alpharr = dindgen(401)/100.d ; Grid from 1 to 5 in 0.01 steps
  ENDIF ELSE BEGIN
     xlim = 0b
     xmax = !values.d_infinity
  ENDELSE
     
  IF (n_elements(vec) EQ 1)  || $
     (n_elements(vec) GT 1 && min(vec) LE 1) THEN BEGIN
     message,'Warning: RANGE input must must contain a vector; '+$
             'using default.',/inf
     undefine,vec
  ENDIF
  
  IF (n_elements(sample) GT 1) || $
     (n_elements(sample) EQ 1 && sample LT 2) THEN BEGIN
     message,'Warning: SAMPLE input must be a positive integer > 1; '+$
             'using default.',/inf
     undefine,sample
  ENDIF
  
  IF (n_elements(limit) GT 1) || $
     (n_elements(limit) EQ 1 && limit LT min(x)) THEN BEGIN
     message,'Warning LIMIT input must be a positive value >= 1; '+$
             'using default',/inf
     undefine,limit
  ENDIF
  
  IF (n_elements(xminx) GT 1) || $
     (n_elements(xminx) EQ 1 && xminx GE max(x)) THEN BEGIN
     message,'Warning: XMIN input must be a positive value < max(x)'+$
             'using default.',/inf
     undefine,xminx
  ENDIF
  
  ;; Select Method (discrete or continuous) for Fitting
  ;;   Determine if elements in X are integers
  IF max( abs(x - round(x) ) ) LE 1.d-3 THEN f_dattype = 'INTS' ELSE $
     IF size(x,/type) EQ 4 || size(x,/type) EQ 5 THEN f_dattype = 'REAL' ELSE $
        f_dattype = 'UNKN'
  IF strcmp(f_dattype,'INTS') && min(x) GT 1000 && n_elements(x) GT 100 THEN $
     f_dattype = 'REAL'
  
  f_dattype = 'REAL'            ; Just set to 'REAL'
  ;;=====================================================================
  ;; Estimate xmin and alpha accordingly
  CASE f_dattype OF
     
     'REAL': BEGIN
        
        IF logdel THEN BEGIN    ; Set XMINS by Delta(Log X) scheme
           xs = x[uniq(x,sort(x))]
           xr = alog(minmax(xs[0:n_elements(xs)-2]))
           xmins = exp(dindgen(ldmax+1)/double(ldmax)*(xr[1] - xr[0]) + xr[0])
        ENDIF ELSE BEGIN        ; Set XMINS by original method
           xmins = x[uniq(x,sort(x))]
           xmins = xmins[0:n_elements(xmins)-2]
        ENDELSE
        
        IF n_elements(xminx) THEN $
           xmins = xmins[(where(xmins GE xminx))[0]]
        IF n_elements(limit) THEN $
           xmins = xmins[where(xmins LE limit)]
        IF n_elements(sample) THEN $
           xmins = xmins[uniq(round(linspace(0,n_elements(xmins)-1,sample)))]
        
        dat = dblarr(n_elements(xmins))
        z   = x[sort(x)]
        FOR xm=0,n_elements(xmins)-1 DO BEGIN
           xmin = xmins[xm]
           z    = z[where(z GE xmin)]
           n    = double(n_elements(z))
           
           ;; Estimate alpha using direct MLE
           IF xlim THEN BEGIN
              xr    = xmin/xmax
              xra   = xr^(alpharr - 1.d)
              dLda  = n / (alpharr - 1.d) - total( alog(z/xmin), /DOUBLE ) + $
                      n * alog(xr) * xra / (1.d - xra)
              dLda2 = n / (alpharr - 1.d) - total( alog(z/xmin) )
              sind = sort(dLda)
              IF oldver THEN $
                 a = interpol(alpharr[sind],dLda[sind],0.d) - 1.d ELSE $
                    a = interpol(alpharr[sind],dLda[sind],0.d, /NaN) - 1.d
              ;; cgPlot,alpharr,dLda,yr=[-1,1]*20,tit=strtrim(xmin,2),thick=2,$
              ;;        psym=-16
              ;; cgOplot,alpharr,dLda2,color='orchid',thick=2,psym=-16
              ;; vline,a+1.d
              ;; vline,0.,/h,color='blue'
              ;; vline,n / total( alog(z/xmin) ) + 1.d,color='orchid'
              ;; wait,0.05
           ENDIF ELSE a = n / total( alog(z/xmin) )
           
           IF nosmall THEN BEGIN
              IF (a-1)/sqrt(n) GT 0.1 THEN BEGIN
                 dat = dat[0:xm-1]
                 xm = n_elements(xmins)+1
                 BREAK
              ENDIF
           ENDIF
           
           ;; Compute KS statistic
           cx      = dindgen(n)/n
           cf      = 1.d - (xmin/z)^a
           dat[xm] = max( abs(cf-cx) ,mazi)
           ;; cgPlot,z,cx,thick=2,/xlog,psym=10
           ;; cgOplot,z,cf,color='orange',thick=2,psym=10
           ;; vline,z[mazi],thick=2,color='teal'
           ;; print,xmin,xmax,xr,a,dat[xm]
           ;; wait,0.2
        ENDFOR                  ; End loop over XMINS
        
        save,xmins,dat,filename='./masses/save_files/ksd.sav',VERBOSE=verbose
        ngt = reverse( findgen(n_elements(xmins)) + 1.)
        dat2 = dat / (1.22/sqrt(ngt))
        D     = min(dat2,D_min)
        xmin  = xmins[D_min]
        z     = x[where(x GE xmin)]
        n     = double(n_elements(z))
        IF xlim THEN BEGIN
           xr = xmin/xmax
           dLda = n / (alpharr - 1.d) - total( alog(z/xmin) ) + $
                  n * alog(xr) * xr^(alpharr-1.d) / (1 - xr^(alpharr-1.d))
           sind = sort(dLda)
           IF oldver THEN $
              alpha = interpol(alpharr[sind],dLda[sind],0.d) ELSE $
                 alpha = interpol(alpharr[sind],dLda[sind],0.d,/NaN)
        ENDIF ELSE alpha = 1.d + n / total( alog(z/xmin) )
        IF finite THEN alpha = alpha * (n-1)/n + 1/n ; Finite-size correction
        IF n LT 50 && ~finite && ~nowarn THEN $
           message,'Warning: finite-size bias may be present.',/inf
        L = n*alog((alpha-1)/xmin) - alpha * total( alog(z/xmin) )
        
        ;;print,n,finite
     END
     
     
  ENDCASE
  
  IF xlim THEN x = retx
  
  RETURN, alpha
END
